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Abstract. A new numerical algorithm for calculating isother- 
mal wind accretion flows has been developed and is applied 
here to the analysis of the hydrodynamics of two-dimensional 
plane symmetric accretion flows in wind-fed sources. Polar co- 
ordinates are used to ensure fine resolution near the object. It 
is found that a thin accretion column is formed which shows 
wave-like oscillations. Small accretion disks are formed tem- 
porarily around the object. Mass accretion rate and angular mo- 
mentum accretion rate exhibit quasi-periodic oscillations. The 
amplitudes of the oscillations depend on the size of the inner 
boundary, the number of grid points and the method of calcu- 
lation. For a smaller size of the inner boundary, finer grids and 
more accurate numerical schemes, the amplitudes of the oscil- 
lation become larger. 

Key words: Accretion, accretion disks - Hydrodynamics - In- 
stabilities - Shock waves - Methods: numerical - pulsars: gen- 
eral 



1. Introduction 

Wind-fed accretion by a compact gravitating object is an im- 
portant astrophysical phenomenon, which arises in many situ- 
ations including massive X-ray binaries. 

Wind accretion has been studied exte nsively both analyti- 
cally and numerically. Hoyle & Lyttleton ( 1939[ ) were the first 
to study an axisymmetric wind accretion onto a gravitating ob- 
ject immersed in a uniform flow. Assuming that the fluid par- 
ticles collide at the axis and lose their tangential momentum, 
they obtained a mass accretion rate of, 



M = TrpocVocRl 

where the accretion radius, R a , is given by 

2GM 

Ra — 



V 2 



(1) 



(2) 



and poo, Voo, G, M are the density, velocity of the gas at in- 
finity upstream, the gravitational constant and the mass of the 
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object, respectively. In this estimate the pressure is neglected 
which means that the Mach number of the gas is infinitely large. 



Bondi & Hoyle (1944) have worked on the model in more 



detail and Bondi (1952) obtained a formula for the case of 



spherical accretion. The latter case means that the Mach num- 
ber at infinity is zero. He also proposed an interpolation for- 
mula for flows with finite Mach number. 

The numerical computation of axisymmetric wind accre- 



tion was first attempted by Hunt (1971, 1979) and general 



agreement between the accretion rates from analytical esti- 
mates and his numerical calculations was found. Shima et al. 



( |1985| ) and Koide et al. ( Jl 99 1| ) calculated the same problem 
using a finer mesh, a modern algorithm and a supercomputer, 
and found accretion rates about a factor two larger than those 
obtained by Bondi. 

If the flow at large distances is not uniform, angular mo- 
mentum can be transferre d to th e compact object causing spin 
up or spin down. Ho et al. (1989) considered this problem based 
on non-axisymmetric simulations. Two-dimensional (planar) 
numerical computations of an inhomogeneous flow were car- 
ried out and non-steady 'flip-flop' motions were found by sev- 
eral auth ors (M atsuda et al. |1987| , Taam & Fryxell[l988j Fryxell 
& Taam 19 8 8b 



Furthermore it was found that non-steady motion occurs 
even in a homogeneous medium (Matsuda et al. 1991). Livio 



et al. (1991 ) suggested a possible cause of the instability. Bof- 
fin & Anzer (1994 ) calculated two-dimensional wind accretion 



using a smoothed particle method (SPH) and confirmed the un- 
steadiness. Benensohn et al. (1997) calculated uniform wind 
accretion of an adiabatic gas with polytropic index 7 = § and 
got 'flip-flop' nonsteady motions. 

Three-dimensional simulations of wind accretion were per- 
formed and compared with two-dimensional cases by Sawada 



et al. (1989) and Matsuda et al. ( 1991 1992). The three-dimen- 



sional calculations performed by Ishii et al. (1993), Ruffert 
( p94a[ [l994b| , [199^ , [l99§) and Ruffert & Arnett ( |1994| ) have 
sufficient spatial resolution to lead to non-steady oscillatory 
flows. Generally these oscillations have lower amplitudes than 
the corresponding two dimensional flows. 
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Since the real size of any astrophysical object is usually 
many orders of magnitude smaller than the accretion radius, 
it seems impossible to calculate the flow all the way to the 
surface. Thus the inner boundary of the computations is much 
larger than the real object. This is an artifact of the simulations. 

But one finds that most of these calculations show non- 
steady behaviours of various degrees. The flows for which 7 
is close to units tend to exhibit the largest fluctuations. 

In this present investigation we want to focus on the spe- 
cial case of isothermal flows. Wind ac cretion with 7 close to 
unity h as been studied by Matsuda et al ( 1991 ), Boffin & Anzer 
(1994), and Ruffert (1996). Since these flows were more vi- 
olent than those of higher 7 they seem of particular interest. 
Physically they can be thought of as approximations to config- 
uration of low optical depth and efficient radiative losses result- 
ing in a temperature distribution which is almost uniform. 

In principle it would be desirable to include an energy equa- 
tion which describes all gain and loss mechanisms of the flow. 
But this approach would complicate our calculations enor- 
mously. Therefore we shall use the approximation of isother- 
mality. 

The objective of this study is to calculate the wind accre- 
tion flows of an isothermal gas at various Mach numbers, M, 
and to study the unsteadiness of the accretion column and the 
appearance of accretion disks around compact objects using a 
fine mesh and high resolution scheme. 

Here we present the results for two dimensional wind ac- 
cretion of an isothermal gas by a new high resolution numerical 
scheme with fine numerical meshes. 

2. Governing Equations and Numerical Procedure 

2.1. Governing Equations 

The governing equations are the hydrodynamic equations of 
an isothermal gas with gravity. The selfgravity of the gas is 
neglected. These equations are written in integral form as, 

<9Q 



dt 



■dr- 
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, r = \f x 2 + y 2 , 



(3) 



(4) 



(5) 



(6) 



(7) 



where p,u,v,p represent density, velocity in x-y direction and 
pressure, dr and ds are the differential volume and surface el- 
ements and (x n , Vn) is the outward normal of the surface ele- 
ment. The pressure for an isothermal gas is defined by. 



where c is the constant sound speed. 

The two momentum equations can be also written in terms 
of the angular momentum and of the momentum in the local 
radial direction, leading to 



+ f [y(mu + px n ) - x(mv + py n )] ds — 0, (9) 
Jon 



f d 

J — [p{uX r + vY r )} dr 



+ d> [X r (mu + px n ) + Y r (mv + py n )} ds 

J dQ. 
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(11) 



These momentum equations are used in the angular mo- 
mentum conserving (AMC) scheme which will be described 
later. 



2.2. Initial and Boundary Conditions 

Initially the whole computational field is filled with a uniform 
gas. The outer boundary condition is given by the analytic solu- 
tion of Bisnovatyi-Kogan et al. ( 1979[ ). An absorbing boundary 
condition, in which the density is very low and the flow velocity 
is zero, is used at the inner boundary. 



2.3. Computational Mesh 

A two dimensional cylindrical grid is used for our computa- 
tions. The mesh is divided uniformly in angular zones and is di- 
vided nonuniformly in radius. The mesh size in the radial direc- 
tion (Ai?) changes exponentially such as, ARi = Ai?ia l_1 . 

The radius of the absorbing inner boundary (R m i n ) is 0.01 
i? a and that of the outer boundary (i? max ) is 10 R a . The stan- 
dard grid is made of 200 mesh points in the angular direction 
and 140 mesh points in the radial direction. For the fine grid, 
these numbers are 400 and 250, respectively. The finest mesh 
size in the radial direction (Ai? x ) is 0.001 i? a in both the stan- 
dard and the fine grid (see Table |l]). 

Under astrophysical conditions, the inner boundary is usu- 
ally located at the magnetopause of the compact accreting ob- 
ject. This radius is in general much smaller than the inner 
boundary of our computational domain. But even if the two val- 
ues are comparable one has to realize that the numerical bound- 
ary conditions do not describe the physical situation at the mag- 
netopause which will be much more complex. The absorbing 
boundary condition at the inner radius is more or less an ar- 
tificial construction. In order to study the influence of i? m ; n , 



p = pc 



(8) computations with different values of i? m i n were carried out. 
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Fig. 1. Density contour for the AM04CL005 case, at times T=22.1, 22.3, 22.5, 22.7 and 22.9. The scale is logarithmic from to 
3, in steps of 0.3 . 
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Radial Mesh 
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A7?i/i? a 


a 


Standard grid 


200 


140 


10 


0.01 


0.001 


1.04533 


Fine grid 


400 


250 


10 


0.01 


0.001 


1.02200 



Table 1. Computational grids used in the simulations 



2.4. Finite Volume Method 

By applying the hydrodynamic equations in integral form to a 
rectangular volume constructed by mesh lines, a finite volume 
formulation is obtained. In order to obtain high spatial resolu- 
tion, a MUSCL type approach and van Albada's flux limiter 
are used. The numerical fluxes which are used by the MUSCL 
scheme are described in section 2.5. 

Hydrodynamic and gravitational terms are integrated si- 
multaneously. A two-step Runge-Kutta method is used for time 
integration. This scheme is accurate to second order in both 
space and time. 



2.5. Simplified Flux Splitting ( SFS) for Isothermal Gas 

The finite volume method for solving the hydrodynamic equa- 
tions takes the approximation of the flux normal to the interface 
which are denoted by Eqs.( |J),( ||) and ( |^). 

The numerical flux in the MUSCL scheme is computed 
from the values defined on both sides of the cell interface by, 

F = F(Q+,Q_) (12) 

where ± denote the values on the left(+) and right(— ) side, and 
m is the mass flux. This flux is computed by solving the Rie- 
mann initial value problem using the left and right side physical 
values. It is in principle possible to solve this problem exactly, 
but this would require a very large amount of computer time. 
Therefore approximate but very fast algorithms have been de- 
veloped. 
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Fig. 2. Density contour for the AM04CL005 case, at times T=22., 22.2, 22.4, 22.6 and 22.8. The scale is logarithmic from to 3, 
in steps of 0.6, and we have zoomed close to the star. The arrows visualize the velocity field 



The way in which these are calculated is vital for the 
MUSCL type scheme, since it influences the computer time, 
robustness and accuracy. Both the robustness for strong shock 
waves and the accuracy for slip surfaces are necessary for this 
study, because of the high mach number near the object and the 
strong shearing in the accretion column and the disk. 

Many algorithms for the numerical flux have been devel- 
oped. Approximate Riemann values of the fluxes were ob- 
tained by the Roe ( ]198l[ ) and (Chakravarthy & Osher |l982j ). 
Although they are exact for shocks or rarefactions of moderate 
strength, they give unphysical results for strong rarefactions or 
shocks respectively. Thus they are not robust enough for hy- 
personic computations like this study where very strong shocks 
and rarefactions are produced by the high speed flow acceler- 
ated by gravity. Flux vector splitting schemes (Steger & Warm- 



the supersonic accretion flow. Wada & Liou (1994) discussed 



ing 1981 van Leer 1982 Hanel & Schwane 1989 ) are simple 



and more robust, and it has been shown that they have enough 
accuracy for the shock tube problem or flows around airfoils. 
However they are not accurate for contact discontinuities such 
as slip surfaces, since they have excess numerical shear stress 
which influences the momentum transfer especially in accre- 
tion disks. Thus these existing schemes are not good enough for 



these problems in the context of to aerospace applications. Note 
that we consider only discontinuities in one-dimension or those 
aligned to grid lines in multidimensional flows. In general, dis- 
continuities do not align to grid lines and all schemes including 
exact Riemann solver can not capture them exactly. Neverthe- 
less it has been recognized that those characteristics strongly 
affect the accuracy. 



Recently Liou & Steffen (1993) developed an Advection 
Upstream Splitting Method (AUSM). Althogh the AUSM ex- 
hibits small overshoots at shocks, it is as simple and robust 
as flux vector splitting schemes and as accurate as the exact 
solver for contact discontinuities. Inspired by this work, a fam- 
ily of AUSM type schemes have been developed (Jounouchi 
et al. Il993l Wada & Liou [19941 Shima & Jounouchi Il997|). It 



has been shown that these AUSM type schemes are sufficiently 
simple robust and accurate. 



Jounouchi et al. (1993, see also Shima & Jounouchi 1997) 



developed a new AUSM type numerical flux scheme named 
Simplified Flux Splitting (SFS). They reduced the overshoot at 
the shock wave. Originally the SFS was developed for adia- 
batic gas, but it has been modified for an isothermal gas in this 
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study. The SFS scheme for an isothermal gas can be written as 
follows. 



R=0.005 



Rmin=0.005 



SFS 



m + \m\ m — \m\ 



* +pN, 



(13) 



where p is an average value of the pressure defined by the 
following relations: 



P = P+P+ +/3-P-, 

P± = i(2TM±)(Af±±l) 2 , 

M± = min(l, max(0, M±)), 
and 



(14) 
(15) 
(16) 

(17) 



The mass flux m is given by the flux vector splitting method 
as follows: 

777 = 777 + + 777_, (18) 

m± = ±0.25 p±c{M ± ± l) 2 for\M±\ < 1 



m± = 0.5p±(V n± ±\V n± \) for\M±\>l. 



(19) 



Note that the usage of the flux vector splitting does not degrade 
the accuracy of a contact discontinuity for isothermal gas, be- 
cause at the slip surface the SFS scheme gives zero flux normal 
to the surface. 

2.6. Angular momentum conserving scheme 

Finite volume methods are usually based on the conservation 
laws of mass, linear momentum and total energy. Conserva- 
tion of energy is not used in the present computations because 
the gas is isothermal. These conservation laws, of course, also 
imply the conservation of angular momentum for the exact so- 
lutions. The conservation of linear momentum, however, is just 
an approximation for the angular momentum in a discretized 
description. As the accuracy is at most second order in mesh 
size with any finite volume method, it depends on the numeri- 
cal mesh. 

If the accreting matter in the accretion column does not col- 
lide with itself and if it keeps some angular momentum , it will 
follow a Keplerian orbit and cannot accrete. Thus the forma- 
tion of an accretion disk near the object directly reflects the 
angular momentum of the accreting matter, and it therefore de- 
pends crucially on the accuracy of the conservation of angular 
momentum. 

The angular momentum conserving (AMC) scheme based 
on Eqs. (^J) and ( |Io[ ) shows less mesh dependency than the lin- 
ear momentum conserving (LMC) scheme as will be shown 
later. The AMC scheme can easily be obtained from the LMC 
scheme. In fact, the difference between the AMC and LMC 
computer codes is only a few lines. 
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Fig. 3. Time history of the mass accretion rate (left) and angular 
momentum accretion rate (right) for the AM04CL005 case at 
times corresponding to Fig. 0, AMC scheme 



2. 7. Accurate local time stepping 

The solutions of wind accretion are nonsteady in most cases. 
Therefore, the numerical scheme has to be accurate in time to 
capture the motion correctly. Two-step Runge-Kutta time step- 
ping with a small Courant number (less than 0.2) was used for 
this purpose. Near the object, the flow has very high velocity 
and a fine mesh is required there for spatial resolution, thus, the 
size of the time step is very restricted by the Courant condition 
and very long computing times are required to reach physically 
meaningful solutions, when the usual uniform step is taken. 

In order to avoid this problem, time accurate local time 
stepping was used in this investigation. First a global time step 
is defined, then the time step at each cell is divided into two 
if the Courant number is larger than our limit. This procedure 
is repeated until Courant condition is satisfied in all cells. The 
inner iteration is carried out in one global iteration when the 
local time step is smaller than the global time step The conser- 
vation laws are fulfilled by summing up the fractional fluxes 
correctly. A speed up by about a factor 20 was obtained by 
this procedure. However, over 100 hours of computer time on 
a 2 GFLOPS class supercomputer (Fujitsu VPP-300) were still 
needed for the fine mesh case. 

2.8. Units 

The obvious reference length is the accretion radius i? a . The 
fluid velocity far upstream, Voa , is used as a reference velocity, 
thus the time unit is 



T, 



ref 



(20) 



The mass accretion rate is normalized by the two-dimensional 
Hoyle & Lyttleton rate given by 

Mm = 2p oc V oc R a . (21) 
The angular momentum accretion rate is normalized by 

^rcf = Poo^oo 2 -Ra- (22) 
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Table 2. Parameters and nondimensional values of each case. M 1 is the root mean square of fluctuations from the averaged value 
of M. AVE(M) and AVE(J) are the average of mass and momentum accretion rate respectively, and RMS(J) is the root 
mean square of J. 



3. Numerical Results 

3.1. Calculated Cases 

The cases which we have calculated are summarized in table 
^. Our main interest is in the thin accretion column of super- 
sonic flows, thus Mach number A4=4 was chosen as the stan- 
dard case. The results obtained with the LMC scheme exhibited 
more mesh dependency than those of the AMC scheme, there- 
fore we show mainly the AMC results 

3.2. Mesh Dependency of the LMC and AMC schemes 

Table || summarizes how for both the LMC and AMC scheme 
the mass accretion rate, M , and the RMS value of the angu- 
lar momentum accretion rate, RMS(J) depend on the details 
of the models. The mass accretion rate changes from 0.420 to 
0.817 when the size of the central hole increases from 0.005 
i? a to 0.02 i? a in the LMC calculations (case LM040.005, 
LM040_01, LM040D2). On the other hand, the mass accre- 
tion rates for the AMC cases are in the range 0.838 to 0.951. 

For the LMC scheme, as the central hole becomes smaller, 
mass and angular momentum accretion also become smaller 
with the standard grid. When the central hole is large, a tempo- 
rary accretion disk is formed and the averaged mass accretion 
rate is close to the Hoyle & Lyttleton value. But for a small cen- 
tral hole an almost permanent accretion disk is formed near the 
hole which blocks further accretion. However in the fine grid 
case this decrease of the mass accretion is not found and the 
existence of the accretion disk is temporary. On the other hand 
for the AMC scheme, accretion disks are always transient and 
accretion rates are similar. 

This fact indicates that the conservation of angular momen- 
tum is very important for the accretion problem. Since the re- 
sults based on the AMC scheme have less mesh dependency, 



Rmin=0.02 .„ Rmin=0.01 




50 100 150 50 100 150 



TIME TIME 

Fig. 4. Time history of the mass accretion rate. Al=4.0; 
i? min =0.02, .01, .005 and .01 fine grid, AMC scheme 



they seem to be more reliable. Thus only the results from the 
AMC scheme will be shown in the rest of this section. 
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Fig. 5. Time history of the angular momentum accretion rate. 
M=4.0; R mia =0.02, .01, .005 and .01 fine grid, AMC scheme 



3.3. Occurrence of Oscillations and Formation of Accretion 
Disks 

The accretion column of the supersonic accretion flow is nar- 
row and large oscillations are found. In Figure [j], the typi- 
cal sequence of the swinging of the accretion column of case 
AM040005 is shown. The oscillating accretion column swings 
over 180 degrees. This oscillation is similar to that shown in 
the computation of Boffin and Anzer ( 1994 1, but there is also 
an accretion disk in our computation. The radius of the accre- 
tion disk is about 0.1 i? a . 

Enlarged views are shown in Fig. ||. The formation and de- 
struction of the accretion disk is clearly seen in these figures. 
An anti-clockwise rotating accretion disk is formed at T=22. 
The accreting matter falls from the upper-left direction through 
the accretion column, thus this matter has anti-clockwise an- 
gular momentum and it accelerates the disk. Then the column 
is pushed backward by accreting matter from upstream. When 
the column moves behind the object, the accreting matter has 
clockwise angular momentum (T=22.2). This matter collides 
with the disk and destroys the disk (T=22.4). Then a clockwise 
rotating disk is formed (T=22.6 and following). As shown in 
Fig. [3j mass and angular momentum is accreted when the disk 
collapses. 

The averaged mass accretion rate of this case is 0.855, but 
the instantaneous rate is over 10 times larger than the Hoyle & 
Lyttleton estimate. The averaged momentum accretion is al- 
most zero. However, a fair amount of instantaneous angular 
momentum accretion with altering direction is also found. This 
can be understood as follows: the symmetric accretion column 



Fig. 6. Time history of the mass accretion rate (upper frames) 
and angular momentum accretion rate (lower frames). 7W=4.0 
Rmin = 0.01. The results from the standard grid (left) and from 
the fine grid (right) are shown, LMC scheme 



cannot have angular momentum, but the portions of an asym- 
metric accretion column can have angular momentum. When 
the positive angular momentum portion falls downward, a pos- 
itive accretion disk is formed. This accretion disk loses its an- 
gular momentum and accretes onto the object, when it collides 
with the following portion which has negative angular momen- 
tum. 



3.4. History of Mass Accretion Rate and Angular Momentum 
Accretion 

The histories of mass accretion and angular momentum accre- 
tion of representative cases from the AMC scheme are shown 
in Fig. |] and [|. Those from the LMC scheme are shown in 
Fig. |. 

3.4.1. Effects of i? m i n and Mesh size 

The histories of accretion rates for .M=4.0 from the AMC 
scheme with different i2 m in and mesh size are shown in Figs. |] 
and U As seen in those figures, the fluctuations of the accretion 
rate become bigger with smaller i? m i n . 

We can see from table 2 and Fig. 4 that the mean values 
of the mass accretion rate are quite similar for different values 
of Rmin, whereas the peak amplitudes of Al increase with de- 
creasing Rmin- This then implies that the accretion bursts are 
correspondingly shorter for smaller values of R m in- Although 
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Fig. 8. Time history of the mass accretion rate (left) and of the 
angular momentum accretion rate (right). .M=8.0 i? m in=o.oi 



tions. The conservation of angular momentum is only achieved 
to second order accuraey in the LMC scheme. Thus, we believe 
that the standard grid is not fine enough for the LMC scheme 
to capture the accretion disk formation. 

3.4.2. Effects of the Mach Number 

If the Mach number of the uniform flow is 1, the flow is per- 
fectly steady, but for larger Mach number it is always non- 
steady. Note that the case of A4=IA is also non-steady and 
shows small fluctuations from the beginning. However, non- 
steadiness with large amplitudes occurs only after 40 time 
units. The Mach number dependence of the averaged mass ac- 
cretion rate, the fluctuation of the mass accretion rate and the 
fluctuation of the angular momentum accretion rate are sum- 
marized in Fig. 0. These are the results from the AMC scheme 
with the standard grid. 

As shown in Fig. the amplitude of the oscillations is 
largest when the Mach number is 2, and the fluctuations be- 
come smaller as the Mach number becomes greater than 4. This 
is also shown in the history of the mass and angular momentum 
accretion of the Mach 8 case (case AM08CL01) in Fig. ||. 

As the Mach number increases, the amplitudes of the ac- 
cretion column oscillation become smaller. This behavior is 
clearly shown in the time averaged density contours of Fig. ^. 
In the time averaged solution, the wiggles of the accretion col- 
umn are smeared out and the high density regions around the 
accreting object are seen. The high density regions are concen- 
trated in a narrow cone in the high Mach number cases (see the 
.M=8 case in Fig. [)]). It is interesting that the transient accre- 
tion disks are also seen as high density regions in the averaged 
solutions and that averaged accretion columns look like bow 
shocks. 



the results for the fine grid (Fig. ^] & |J) have fluctuations of 
higher frequency and are more chaotic, they are generally sim- 
ilar to those for the standard grid. The similarity is also clearly 
seen in Table ||. All values except M' are similar for different 
i? m i n . The sequences of the oscillations are also similar. The 
difference in M' can be explained by the mechanism of the 
accretion disk destruction. The orbits of the anti-clockwise ro- 
tating accretion disk in Fig. || are almost circular, but become 
elliptical when the accreting mater which has clockwise angu- 
lar momentum interacts with matter of the anti-clockwise ro- 
tating disk (T=22.2 in Fig. |[). The matter is taken away from 
the computational space when its elliptic orbit touches the cen- 
tral hole. Thus, the bigger central hole absorbs the matter over 
a longer period than the smaller hole. Therefore the fluctuation 
of mass accretion becomes larger for a smaller hole. 

On the other hand, the LMC scheme shows a different i? m i n 
dependence. For the standard grid, the mass accretion rate be- 
comes smaller as i? m i n decreases, as can be seen from the av- 
erage values of M which are given in table 2.. Because an al- 
most permanent accretion disk is formed it blocks further ac- 
cretion. But this permanent accretion disk is not found for the 
fine grid. In this case the accretion rates show large fluctua- 



4. Summary and discussion 

Two dimensional planar wind accretion flows of a supersonic 
isothermal gas were analyzed numerically. The computational 
efficiency was improved as compared to earlier calculations by 
using an accurate local time stepping and a new type of upwind 
scheme. 

Comparing our new results with earlier calculations one 
realises that 2D flows with 7 close to unity have so far not 
been studied very systematically. Matsuda et al. (1992) pre- 
sented one case with 7 = 1.05 and one with 7 = 1.005. Their 
resolution near the central object was much coarser than our 
resolution. They had also obtained flip-flop configurations, but 
the instability took much longer to develop. Their amplitudes 
were smaller and the widths of the accretion columns consider- 
ably larger. Boffin and Anzer (1994) also calculated one model 
with 7 = 1.1 which has a narrow oscillating accretion col- 
umn, but again a low amplitude. 3D flows with 7 = 1.01 have 
been systematically studied by Ruffert ( 1996). In these flows 
the accretion regions are always much wider and do not show 
the systematic oscillations of entire accretion region, although 
the accretion itself is non-steady. There differences between 3D 
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Fig. 9. Time averaged density contour. Rmin=0.01, A4=2,4 and 8. The contours are logarithmically spaced between and 3, in 
steps of 0.3 



and 2D flows are quite general and independent of the value of 
7- 

We also have developed a numerical scheme which ex- 
actly conserves angular momentum (AMC scheme) and com- 
pared its performance to the earlier linear momentum conserv- 
ing (LMC) scheme. We found that the results obtained with 
our LMC scheme depend very strongly on the resolution of the 
grid used, whereas the AMC scheme gives almost no differ- 
ence for the two types of resolution considered here. Therefore 
we feel that the AMC calculations are much more reliable and 
all the results presented in the previous section are based on 
this scheme. The fact that our AMC and LMC results differ 
by large amounts suggests to us that many of the earlier calcu- 
lations which basically conserve linear momentum should be 
taken with caution. This will be particularly important for the 
temporal behavior of the accretion of angular momentum J (t). 

Our calculations of supersonic flows show that in all cases 
the accretion of both mass and angular momentum is very er- 
ratic. It occurs in very short intervals and at extremely high 
rates. This behavior can be explained by the formation of Ke- 
plerian disks near the inner boundary. If the specific angular 
momentum of the infalling material is larger than that of the 
Kepler orbit at the innermost radius then this material cannot 
accrete. 

If such high angular momentum material is flowing in long 
enough, a disk will form which blocks further accretion very 
efficiently. Material with very low angular momentum or with 
opposite rotation falling in during a subsequent phase interacts 
with the disk and can destroy it. This will lead to a burst of 
the accreted mass and angular momentum. After such a burst 
the process can be repeated and a new disk will form. Our cal- 
culations indicate that reversals of the disk rotation are quite 
common. 

For the modeling of X-ray binaries fed by wind accretion 
the fluctuations of J are of major importance. They can be - 
brought into relation with the observed spin- up an d spin-down 
of these X-ray pulsars; see Anzer & Borner ( 1995 ). In their in- 



vestigation they showed that the random fluctuations calculated 
by Ruffert (1994c ) for 3d models were by a factor ps 10 too low 
in order to explain the pulse period variations observed in the 
source Vela X-l. However our new calculations give variations 
of J which are substantially larger than those found by Ruffert. 
We have obtained typical values of RMS (J) of the order of 0.2 
in our dimensionless units (see Table |2|). This result can also be 
formulated as: 



RMS(J) = 0.2 Pc 



Ri 



= 0.1 M HL voo R a 

RJ 0.1 Woo -Ra 

since M is ty pically of the order unity. On the other hand 
Ruffert dl996j ) gives RMS(j)=0.01Voo-R a for 7 = 1.01 and 
RMS(j) = 0.03w oo i? a for 7 = 5/3. Taking into account that 

RMS( j)=RMS (Mj) and M » 1 we have RMS ( j) = (0.01- 
0.03)Wooi? a . Therefore our values for the fluctuations area fac- 
tor 3 to 10 larger than those of Ruffert's 3D calculations. There- 
fore on the basis of our calculations one might conclude that 
the observed period fluctuations could in principle, be caused 
by random fluctuations of J. But the amplitudes are only 
marginally large enough and any slight reduction of the effi- 
cency would rule out this interpretation. There is in particular 
the aspect that our calculations are two-dimensional whereas 
the real flows are three-dimensional and the difference in am- 
plitudes between 2D and 3D flows could be sufficiently large 
to make the described interpretation invalid. To really answer 
this question requires full 3D computations, taking angular mo- 
mentum conservation into account. 
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